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I. INTRODUCTION 

Pairing correlations affect properties of atomic nuclei 
in a profound way 0, 0, 0, 3 • They impact nuclear bind- 
ing, properties of nuclear excitations and decays, and 
dramatically influence the nuclear collective motion. In 
particular, pairing plays a crucial role in exotic, weakly 
bound nuclei in which the magnitude of the chemical po- 
tential is close to that of the pairing gap. In such systems, 
a naive independent single-particle picture breaks down 
and the pair scattering, also involving the continuum part 
of the phase space, can determine the very nuclear exis- 
tence la]. 

Many aspects of nuclear superfluidity can be success- 
fully treated within the independent quasiparticle frame- 
work by applying the Bardecn-Cooper-Schrieffer (BCS) 
m or Hartree-Fock-Bogoliubov (HFB) approximations 
Pf. The advantage of the mean-field approach to the 
pairing problem lies in its simplicity that allows for a 
straightforward interpretation in terms of pairing fields 
and deformations (pairing gaps) associated with the 
spontaneous breaking of the gauge symmetry. However, 
this simplicity comes at a cost. In the intrinsic-system 
description, the gauge angle associated with the particle- 
number operator is fixed; hence, the particle-number in- 
variance is internally broken Therefore, to relate 

to experiment, the particle- number symmetry needs to 
be, in principle, restored. 

Some observables, like masses, radii, or deformations 
are not very strongly affected by the particle-number- 
symmetry restoration, while some other ones, like even- 
odd mass staggering or pair-transfer amplitudes are influ- 
enced significantly. Moreover, quantitative impact of the 
particle-number projection (PNP) depends on whether 
the pairing correlations are strong (open-shell systems) or 
weak (near closed shells). Therefore, methods of restor- 



ing the particle-number symmetry must be implemented 
in studies of pairing correlations. This can be done on 
various levels 2, 7], including the quasiparticle random 
phase approximation, Kamlah expansion [a, la, Lipkin- 
Nogami (LN) method [lS|ll|llEllIll3[l»,']^,'ll, 
the particle-number projection after variation (PNPAV) 
the projected LN method (PLN) [H |Tg|, |20i, Ip 
|22| . and the variation after particle-number projection 

(vAPNP) [Ulil Hill 113. 

In this work, we concentrate on the VAPNP method. 
Recently, it has been shown that the total energy in 
the HFB-I- VAPNP approach can be expressed as a func- 
tional of the unprojected HFB density matrix and pairing 
tensor. Its variation leads to a set of HFB-like equations 
with modified self-consistent fields. The method has been 
illustrated within schematic models [271, and also im- 
plemented in the HFB calculations with the finite-range 
Gogny force 0,115 ■ Here, we adopt it for the Skyrme 
energy-density functionals and zero-range pairing forces; 
in this case the building blocks of the method are the lo- 
cal particle-hole and particle-particle densities and mean 
fields. In the present study, the HFB equations are solved 
by using the Harmonic Oscillator (HO) basis, but the for- 
malism can be straightforwardly applied with the Trans- 
formed Harmonic Oscillator (THO) basis 0,113, which 
helps maintain the correct asymptotic behavior of the 
single-quasiparticle wave functions. 

The paper is organized as follows. Section ITTl erives a 
brief overview of the HFB theory and defines the densities 
and fields entering the formalism. Section llTTl extends the 
VAPNP method of Ref. 21] to the case of the HFB theory 
with Skyrme interaction. The technical details pertain- 
ing to the Skyrme HFB-I- VAPNP method are given in 
Sec. IIVI while Sec. IVl contains an illustrative example of 
calculations for the Ca and Sn isotopes. In particular, 
the LN and PLN approximations are compared to the 
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VAPNP results. Summary and discussion are given in 
Section IVIl Preliminary results of our VAPNP calcula- 
tions were presented in Ref. [29l |. 



II. THE HFB METHOD 

The many-body Hamiltonian of a system of fermions 
is usually expressed in terms of a set of annihilation and 
creation operators (c, c^): 



t t 



nn' mm' ^n^n' m^m^ 



nn'mm' 



where 



Vnn'mm' = (nn'\V\mm' ~ m'm) 



(1) 



(2) 



are the anti-symmetrized two-body interaction matrix- 
elements. 

In the HFB method, the ground-state wave function 
is the quasiparticle vacuum |(f>), defined as Q;fe|<i>) = 0, 
where the quasiparticle operators (a, a^) are connected 
to the original particle operators via the Bogoliubov 
transformation 



* t \ 



Oik = ^ {U*kCn + V, 
n 

^ {VnkCn + Unkci) 



(3) 
(4) 



where the matrices U and V satisfy the unitarity and 
completeness relations: 

U^U + V^V ^ I, UU'' + V*V'^ ^ I, (5) 
U^V + V^U ^0, UV'' + V*U^ ^0. (6) 

A. The HFB equations 

In terms of the density matrix p and pairing tensor k, 
defined as 

p^V*V^, K^V*U'^^-UV\ (7) 
the HFB energy is expressed as an energy functional: 



E[p, k] 



($1$) 



- Tr[(e + ir)p] -iTr[A«:*], (8) 



where 



n' in' 



' Pm'n' 



nn' mm' '^mrn' 



(9) 
(10) 



The variation of the HFB energy © with respect to p 
and K yields the HFB equations: 



H 



Uk 
Vk 



Uk 



where 



n = 



r- A A 

-A* -(e + r)* + A 



(11) 



(12) 



Uk and Vk are the fcth columns of matrices U and V, 
respectively, and Ek is a positive quasiparticle energy 
eigenvalue. Since the HFB state |$) violates the particle- 
number symmetry, the Fermi energy A is introduced to 
fix the average particle number. 



B. The Skyrme HFB method 

For the zero-range Skyrme forces, the HFB formalism 
can be written directly in the coordinate representation 
[3(1 .31. ,32] by introducing particle and pairing densities 



p{ra,r'(j') = ^p{r,r')5aa' 

+ \Y,{a\aM')P:{r,r'), (13) 

i 

p{ra,r'a') = \p{r,r')5aa' 

+ i^(a|a,|a')AKr,r'), (14) 

i 

which explicitly depend on spin. The use of the pairing 
density /5, 



p(r(T, r'cr') ~ — 2(t' K{r, a, r' , — a'), 



(15) 



instead of the pairing tensor k is convenient when re- 
stricting to time-even quasiparticle states where both p 
and p are hermitian and time-even |3lj | . 

The densities p and p can be expressed in the single- 
particle basis: 

p{ra,r'a') = ^ p„„. V;^, (r', cr')V'n(r, cr), (16) 

nn' 

p{ra,r'a') = J^pnn' K'ir' ,a')Mr,<y), (17) 

nn' 

where p„'„ and /5„/„ are the corresponding density ma- 
trices. In this study, we take ipni^, <j) as a set of the HO 
wave functions. 

The building blocks of the Skyrme HFB method are the 
local densities, namely the particle density p{r), kinetic 
energy density T(r), and spin-current density Jy (r): 

p{r) = p(r,r), 

T(r) = VrS/r'Pir,r')l,^^ , (18) 
^vir) = h iy^-^'i)pAry)\r'=r ' 
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as well as the corresponding pairing densities p{r), f{r) 
and ^ij{r). 

In the coordinate representation, the Skyrme HFB en- 
ergy (IHJ can be written as a functional of the local particle 
and pairing densities: 



(^\H\^) f 



(19) 



The energy density 7i(r) is a sum of the particle H{r) 
and pairing energy density H{r): 



H{r)^H{r)+H{r) 



(20) 



The derivatives of E[p, p\ with respect to density matrices 
p and p define the self-consistent particle (h) and pair- 
ing ih) fields, respectively. The explicit expressions for 
H{r), H(r), h{r,(T,a'), and h{r,a,a') have been given 
in Ref. |3l| and will not be repeated here. 

The Skyrme HFB equations can be written in the ma- 
trix form as: 



h-X h 
h -h + X 



(21) 



where 



dE[p,p] 

dpn'n 



^ y dr rnir,a)h{r,a,a')i^n'{r,a'), (22) 



and 



hrin.' 



dE[p,p] 



dpn'n 

= J2j dr r„{r,a)h{r,a,a')^n'ir,a'), (23) 

crcr' 

and ipi^k and f2,k are the upper and lower components, 
respectively, of the quasiparticle wave function corre- 
sponding to the positive quasiparticle energy Ek- After 
solving the HFB equations (|21|l . one obtains the density 
matrices, 



Pnn' — ^ f2,nkV'2,n'kJ 
Ek>Q 

Pnn' = - ^ </'2,nfc¥'l,„'fe, 
-Efc>0 



(24) 
(25) 



which define the spatial densities (|16|l and H17() . 

We note in passing that the derivation of the 
coordinate-space HFB equations 31] is strictly valid only 
when the time-reversal symmetry is assumed. When 
the time-reversal symmetry is broken, one has to intro- 
duce additional real vector particle densities s, j, T [33l |. 
while the pairing densities acquire imaginary parts; see 
Ref. jSJ] for complete derivations. 



III. VARIATION AFTER PARTICLE-NUMBER 
PROJECTION 

A. The HFB+VAPNP method 

It has been demonstrated that the HFB+VAPNP 
energy, 



E''[p,K] = 



($|pJV|$) 

/rf0($|e»"^(*-^)|$) ' 



(26) 



where P is the particle-number projection operator. 



1 



P- ^ — I a0 e 
Zn 



i<l>{N-N) 



(27) 



can be written as an energy functional of the unprojected 
densities ((TJ. 

The variation of Eq. (jSEl results in the HFB-^VAPNP 
equations: 



N 



N ( \ _ 



= 8k 



N 



where 

n 



k \ yN 



-(A^)* -(e^ + r^ + A~)* 



(28) 



(29) 



Equations (|28|l and (|29|l have the same structure as 
Eqs. Hll|) and H12|l . except that the expressions for the 
VAPNP fields are now different IS El, i.e., 



-JV _ 1 



2 , dcj) 2/(0) {r(0)Tr[ep(0)] 
+ [l-2ie-*'^sin</)p(</>)]eC(</))} -Hh.c, (30) 

^ J fd.p y{c^) {r(<^)Tr[r(0)p(0)] 



+ 2[l-2ie-*'^sin#(0)]r((/))C(0)} -Kh.c.,(31) 
= i / d0 y{^)e-^^^C (4,) A{4,) - {...f, (32) 



A^ = 



with 



-i J dc^ y(0){y(</))Tr[A(</))7«*(0)] 
- Aie-"^ sin 4> C(0) A((/))k* (0) } + h.c. 

rnm(^) ~ ^ ^ ^nn' mm' Pm' n' (^4^) ^ 
ii' m' 

^nn' (^4^^ — 2 ^ ^ ^nn' mm' ^mim' (^4^^ ■} 



(33) 

(34) 
(35) 
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where, using the unit matrix /, 



p{4>) = C{4>)p, (36) 

^(0) = C{(I>)k, (37) 

-K{(j)) = e^''>'C\4>)K, (38) 

C{4>) - e^^* [i + p(e2'*_i)]-\ (39) 
1 e-*'^^det(e*'^/) 



x(0) = 



2vr ^det C(0) 

/ d(t)'x{4>'y 



(40) 
(41) 



and 



r(0) = ie~*'^sin<?!)C(0) 



yy(0')e""^'sin</>' C(0')- (42) 



After solving the HFB+VAPNP equations one ob- 
tains the intrinsic density matrix and pairing tensor: 

p^{V'')*{V''f, (F^)*(C/^)^. (43) 

Finally, the total HFB+VAPNP energy is given by 

E^[p, «] = / d0 y(0) Tr (ep(</>) + ir(</>)p(</>)) 

d^y(0) iTr(A((/.)«;*(0)). (44) 

The quantity y{(j)) plays a role of an A^-dependent metric. 
The integrands in Eqs. (|2ni)-(|S21l take the famihar HFB 
limit at (/!)=0, while the integrand in (|33|) vanishes (A^ 
does not appear in the standard HFB approach). 

B. The Lipkin-Nogami method 

The LN method [13, [ill constitutes an astute and effi- 
cient way of performing an approximate VAPNP calcula- 
tion. It can be considered|3 as a variant of the second- 
order Kamlah expansion [a, [3 ' which the VAPNP en- 
ergy (|26|) is approximated by a simple expression, 



Ei^^ = E[p,p\ ~ \2{{N^) ~ N^), 



(45) 



with A2 depending on the HFB state [$) and representing 
the curvature of the VAPNP energy with respect to the 
particle number. The role of A2 in the Kamlah and LN 
methods differs. In the former, A2 is varied along with 
variations of the HFB state |$), while in the latter, this 
variation is neglected. Had the second-order Kamlah ex- 
pression H45|l been exact, the variation of A2 would have 
been fully justified and the method would be giving the 
exact VAPNP energy. However, since the second-order 
expression is, in practical applications, never exact, it 
is usually more reasonable to adopt the LN philosophy, 
in which one rather strives to find the best estimate of 



the curvature A2 instead of finding it variationally in an 
approximate way. 

When the HFB method is applied to a given Hamilto- 
nian, values of A2 can be estimated by calculating new 
mean-field potentials, F' and A', that are analogous to 
the standard mean fields of Eqs. JHl and (fTn|l : see, e.g., 
Refs. [3,E|. However, apart from studies based on the 
Gogny Hamiltonian ,19,] , such a formula was not used, be- 
cause most often the self-consistent calculations are per- 
formed within the density functional approach or by us- 
ing different interactions in the particle-hole and particle- 
particle channels. Moreover, in most studies, such as 
those of Ref. the terms in A2 originating from the 
particle-hole channel are simply disregarded. 

Similarly, as in our previous study 'l8|, here we adopt 
an efficient phenomenological way of estimating the cur- 
vature A2 from the seniority-pairing expression, 

^ ^ Geff Tr'(l - p)K Ti'pn - 2 Tr(l - pfp^ 
' 4 [Trp(l-p)]'-2Trp2(l-p)2 ' 

where the effective pairing strength, 



Geff 



E, 



pair 



is determined from the HFB pairing energy, 

£^pair = -^TrAK*, 

and the average pairing gap |42| . 



A 



Tr'Ap 
Trp 



(47) 



(48) 



(49) 



Expression (|46|l pertains to a system of particles occupy- 
ing single-particle levels with fixed (non-self-consistent) 
energies and interacting with a seniority pairing interac- 
tion. In our method, this expression is used to probe the 
density of self-consistent energies that determine the cur- 
vature A2. All quantities defining A2 in Eq. depend 
on the self-consistent solution and microscopic interac- 
tion, while the effective pairing strength Gcft is only an 
auxiliary quantity. The quality of the prescription for 
calculating A2 can be tested against the exact VAPNP 
results (see Sec. 01. 



C. The Skyrme HFB+VAPNP method 



Following the VAPNP procedure of Sec. 1111 Al one can 
develop the Skyrme HFB+VAPNP equations by intro- 
ducing the gauge-angle-dependent transition density ma- 
trices: 

p(rcr, r'cr', 0) = ^p„„/(0) V,*'('^':Cr')V'n(T-,cr),(50) 

nn' 

p(rcr,rV',0) = ^/5„„'(0) V;^,(r',cr')V'n(T-,cr).(51) 
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In the above equation, the density matrix p„„/ (</)) is given 
by Eq. ^ while 

= e-^^C(0)p. (52) 

The associated gauge-angle-dependent local densities 
p{r,(j)), T{r,(j)), ^^Jir,(j)), p(r,0), f(r,0), and Jjj(r,0) 
are defined by Eqs. (|18|l in terms of the density matrices 
(|50|l and (|51|l . Using the Wick theorem for matrix ele- 
ments 10), one can show that the gauge-angle-dependent 
transition energy density H(r, (f>) can be obtained from 
the intrinsic energy density H^r) simply by substituting 
particle (pairing) local densities with their gauge-angle- 
dependent counterparts (e.g., p{r)^p{r,(j)j). 

In the case of Skyrme functionals, the HFB+VAPNP 
energy H26|l can be expressed through an integral 



where the transition energy reads 



i?(</>) 



= J dr H{r, 



(53) 



0). (54) 



The projected energy (|53|) is a functional E^[p,p] of the 
matrix elements of intrinsic (i.e., (/>=0) matrices p and p. 

In order to compute the derivatives of {p, p) with 
respect to p and p, one should take first the derivatives 
of E^[p,p] with respect to p((/)) and p{(j)), and then the 
derivatives of p{(j)) and p{(j>) with respect to the intrinsic 
densities p and p. For example, 

dE^[p,p] f 

■ = d(t> y{(l)) 



dpn 



1 dv{4>) 



E 

a/3 

E 



_y{4>) dp„ 

dE{(t>) dpo,p{cj)) 
dpap{4>) dp 

dE{(t>) dpo,0{4>) 
dpapi't') dp 

nn' 

dE{c^) dplpi-c^) 



y 



(55) 



With the use of the identity: 

Smm' - 216^*"^ sin Pmm'i4') 

the partial derivatives in Eq. H55|l can easily be calcu- 
lated: 



'W, (56) 



dpnn' 
dpnn- 

dpmm' {4>) 
dpnn' 

dpmm- {4>) 
dpnn' 



y(0) Ynn' {(!>), 



(57) 



2ie-'* sin((/.)p„,™, (0)C„„ (</)) , (58) 

-2ie-"f' Sm{^)pr,'rn'WCmnW, (59) 



(60) 



where n and s„ (snS^ = 1, = — s„) are defined using 
the time-reversal operator T, as 

fipnir, cr) = Sni^n{r, a). (61) 

By inserting Eqs. (|57(l - (|60|l in Eq. (|55|l . the latter reads 



9p 



= / d(t> y(0) r(0) Ei^) 



where 



d0 2/(0) C(0);i(0)C(0) (62) 

d0 2/('/>) 2ie-*"* sin(0) p{(f>)h{(j))C {(f>) +h.c. 

ftrm' (0) = ;^ 7-T 

Op„'n{(P) 

J d'rrn{r<j)h{r,a,a'A)i;^,{ra'), (63) 

fi-rm' ((/>) = 7-T 

OPn'nW 

= J2 fd^'r rn ira)h{r, a, a', 0)^„- (ra'). (64) 

crcr' 

The derivative of E^{p,p) with respect to p can be 
computed in a similar manner. The 0-dependent fields 
h{r,a,a' , 4>) and h{r,a,a' , (j)) are obtained by substitut- 
ing the local particle and pairing densities in the intrinsic 
fields h{r,a,a') and h{r,a,(j') with their gauge- angle- 
dependent counterparts. 

The Skyrme HFB+VAPNP equations can finally be 
written in the form 



(65) 



with particle-hole and particle-particle Hamiltonians 
= I dci,y(4,)Y{cl>)E{4>) 
dc^y{cf>)e-^'* C{<f)h{4>)C{(j,) (66) 
d(t)y{(t))2ie-"'' sin{(l))p{(l))h{(t))C{(l)) + h.c. , 



d4>y{^)e-^^ \h{4>)C{4>) + (...)' 



(67) 



Finally, solutions of the HFB+VAPNP equations (jHSJ 
allow for calculating the intrinsic density matrices as. 



E JV ,N* 



Ek>0 



Pnn 



E N ,N* 



(68) 
(69) 



Bfc>0 
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Let us re-emphasize that the densities and fields that 
enter the Skyrme HFB+VAPNP equations are immedi- 
ate generahzations of the analogous quantities that ap- 
pear in the standard Skyrme HFB formalism. Of course, 
due to the presence of C{(f)) and integrations over the 
gauge angle, the Skyrme HFB+VAPNP calculations are 
appreciably more involved. 



IV. SKYRME HFB+VAPNP PROCEDURE: 
PRACTICAL DETAILS 



A. Two kinds of nucleons 

As one is dealing with Z protons and N neutrons, two 
gauge angles, (j)n and (f>p, must enter the number projec- 
tion operator: 



1 

2^ 



1 

2^ 



(70) 



Consequently, the total projected energy (|53|l becomes a 
double integral, 

E'^ = j d4>n d4>p yn{4>,i) yp{4>p) E{(j)n,4>p), (71) 

where the transition energy density 

E{<t)n,c^p) ^ J dr H{r,(j)„,q^p) (72) 

depends on both gauge angles 0„, 4>p. 

To simplify notation, we use the isospin label 
Q=T3 (?=+! for neutrons and -1 for protons) and 
q=—q. In the following, we shall employ the conven- 
tion y{(j)q)=yq{(j)q), C q)=C q{(t) q) , aud Y {(t>q)=Yq{<t)q) . 

The isospin-dependent particle-hole and particle-particle 
fields l|nni), can be written as: 

hN _ 



d<Pq y{(l)q) 
Y{^q) y{4>q) E{cf,q,cl)q)dcf>q 

e-2»*' C{<l>q) y{c^q) hq{c^q,cj,q)dcl)q ) C{<l>q) 

d(Pq y{4)q) 2ie~"'''' SUl{(t)q)pq{(t)q) 
y{<Pq)hq{(f>q,(l)q)d(f>q^ C{(j)q) + h.C 

d^q y{c^q) e-*^' 

y{cbq)hqi'f>,^^'i)d<l>q) C{^q) + (...)^ 



(73) 



(74) 



In numerical applications, the two-dimensional integrals 
over the gauge angles are replaced by a sum over L„ x Lp 
points using the Gauss-Chebyshev quadrature method 



B. Canonical representation 

The canonical-basis single-particle wave functions, 

XtJ.{r,(T) ^^Wn^, i^n{r,a), (75) 

n 

are defined by the unitary matrix W which diagonalizes 
the density matrices. 



n' 



(76) 



where are the occupation probabilities < < 1 and 
v'^ + uj^ — 1. In the canonical representation, the gauge- 
angle-dependent matrices become diagonal with the di- 



agonal matrix elements given by: 



CM,) 



Pmi't^q) 



"/jg ^ '^/j.q 
ulq + e'^"''vlq ' 



PMyfll - y2 ,g2»,^,^2 



(77) 
(78) 
(79) 



and the determinant of matrix C{(t)q), needed in Eq. H4U|I. 
becomes a product of the diagonal values H77|) . The use 
of the canonical representation significantly simplifies cal- 
culations of the projected fields. 



C. Intrinsic average particle number in the 
HFB+VAPNP method 

The HFB state |$) is a linear combination of eigen- 
states |A^) of the particle-number operator, i.e.. 



|<i>)=^a^|7V), 



N 



where 



|iV) 



(80) 



(81) 



and N\N) = N\N). The HFB+VAPNP method is based 
on the variation of the projected energy (jSEJ, which is 
the average value of the Hamiltonian on the state |iV), 
~ {N\H\N). Obviously, the projected energy does 
not depend on amplitudes oat, although the intrinsic av- 
erage number of particles. 



N ^ ($|7V|$) 



N 



\aN\^N = Tvp, 



(82) 



does depend on a^- 
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The HFB+VAPNP variational procedure gives, in 
principle, the same value of the projected energy inde- 
pendently of the value of N. This independence can, 
however, be subject to numerical instabilities whenever 
the amplitude a^r, corresponding to the projection on N 
particles, is small. Therefore, for practical reasons, one 
is interested in keeping the average number of particles 
N as close as possible to N, which guarantees that am- 
plitude UN is as large as possible. 

In the standard HFB equations 1)21(1 , the average num- 
ber of particles is kept equal to N by adjusting the 
Lagrange multiplier A. However, in the HFB+VAPNP 
approach, A does not appear in the variational equa- 
tions H65|l . because the variation of the constant term 
A($|7VF^|$)/($|F^|$) = AiV equals zero. Therefore, 
the HFB-t-VAPNP equations (jHS) do not allow for ad- 
justing the average particle number N, which, during 
the iteration procedure, may become vary different from 
N. Moreover, such uncontrolled changes of N from one 
iteration to another may preclude reaching the stable self- 
consistent solution. 

In order to cope with these problems, one can artifi- 
cially reintroduce a constant /i, analogous to the Fermi 
energy A, into the HFB+VAPNP equations (gSJl, i.e.. 



4,, 



(83) 



provided it is equal to zero once the convergence is 
achieved. With this modification, the iterations proceed 
as follows. Suppose that at a given iteration, condition 
N=N is fulfilled. Then, no readjustment of N is neces- 
sary, and in the next iteration one proceeds with /i=0. 

Had such an ideal situation continued till the end, a 
nonzero value of /i would have never appeared, and the 
required solution would have been found. In practice, 
this situation never happens, and at some iteration one 
finds that Trp, i.e., the sum of norms of the second com- 
ponents (p^„j, lEHl) of the HFB+VAPNP wave functions, 
is larger (smaller) than N. In such a case, in the next 
iteration one uses a slightly negative (positive) value of 
/i, which decreases (increases) the norms of <f2nk^ 
decreases (increases) the average particle number in the 
next iteration. Since ^ acts in exactly the same way as 
the Fermi energy does within the standard HFB method, 
the well-established algorithms of readjusting A can be 
used. Moreover, as soon as the iteration procedure starts 
to converge, the non-zero values of ^ cease to be needed, 
and thus n naturally converges to zero, as required. In 
practice, we find that the above algorithm is very useful, 
and it provides the same converged solution with any 
value of iV = iV ± AA^, AA^ being a small integer. 



D. The cut-off procedure for the contact pairing 
force 



When using zero-range pairing forces such as the 
density-dependent delta force, one has to introduce the 



energy cut-off j35|. Within the unprojected HFB cal- 
culations, a pairing cut-off is introduced by using the 
so-called equivalent single-particle spectrum jsj- After 
each iteration, one calculates an equivalent spectrum e„ 
and corresponding pairing gaps A„: 



e„ = (1 - 2P„)^« + A, A„ = 2Eny/Pn{l-Pn), (84) 

where E„ is the quasiparticle energy, A is the chemical 
potential, and P„ denotes the norm of the lower com- 
ponent of the HFB wave function. The energy cut-off is 
practically realized by requesting that the phase space for 
the pair scattering is limited to those quasiparticle states 
for which e„ is less than the cut-off energy ecut (usually 
ecut = 60MeV) 

Obviously, the above procedure cannot be directly ap- 
plied to the HFB+VAPNP method, where the intrin- 
sic quantities, in particular the 'quasiparticle' energies 
, do not have obvious physical meaning. A reason- 
able practical prescription for ecut can be proposed in 
terms of intrinsic {(j) = 0) HFB fields h and h. After each 
iteration of Eq. H65() . the average quasiparticle energies. 



E„ 



yN 



h-\ h 

h -h + X 



yN 



(85) 



together with P„=P^, give the equivalent energies JHU. 
Based on the spectrum of e„, the set of quasiparticle 
states appearing below the cut-off energy can now be 
easily defined. At the same time, the Fermi energy A 
(as an auxiliary quantity) can be recalculated in each 
iteration. 
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FIG. 1: (color online) The neutron equivalent single-particle 
energies JMJ for iV=70 and Z=50 obtained in the HFB+LN 
method (first spectrum) , HFB+VAPNP method using the av- 
erage quasiparticle energies En (second spectrum) 18511 . and 
by using the 'quasiparticle' energies Ej^ calculated in the 
HFB+VAPNP method for different values of the intrinsic neu- 
tron N and proton Z numbers (the remaining five spectra). 
The dashed line indicates the position of the LN Fermi energy 
A. 
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The results of such a procedure are illustrated in Fig.^ 
The left-most spectrum shows the neutron equivalent en- 
ergies obtained within the LN method applied to A^=70 
and Z—50, and the dashed line shows the position of the 
corresponding LN neutron Fermi energy A. For e„ < 0, 
this spectrum is very similar to the HF bound single- 
particle energies of this nucleus. Our method, based 
on the average quasiparticle energies gives almost 

identical negative equivalent energies and quite similar 
positive ones. In particular, for highly positive equiva- 
lent energies, in the region of the cut-off energy ecut = 
60 MeV, similar continuum quasiparticle states appear in 
both methods; this guarantees the correct application of 
the cut-off procedure. The five equivalent spectra shown 
on the right hand side of Fig. ^ were calculated directly 
from the unphysical 'quasiparticle' energies obtained 
for several selected values of the intrinsic particle num- 
bers N and Z. It is obvious that these spectra (even at 
N—7Q and Z—50) bear no resemblance to the real single- 
particle spectra and cannot be used to define the cut-off 
procedure. 



V. SAMPLE RESULTS 

To illustrate the Skyrme HFB+VAPNP procedure, we 
carried out calculations for the complete chain of calcium 
isotopes, from the proton drip line to the neutron drip 
line, and for the chain of tin isotopes with 70 < < 90. 
We used the Sly4 Skyrme force parameterization Isjand 
the mixed delta pairing [ssl l39l | . The calculations were 
performed in the basis of 20 major HO shells. We took 
L=13 gauge-angle points, and this practically ensures ex- 
act projection for all considered nuclei. We have found 
that the HFB+VAPNP procedure is just L-times slower 
compared to the PLN method. 

In our standard HFB calculations js^ls^, the strength 
of the pairing force (assumed identical for protons and 
neutrons) is usually adjusted at a given cut-off energy 
Ecut = 60 MeV to the experimental value of the aver- 
age neutron gap A„=1.245 MeV in -'^^°Sn. In the present 
study, we used this procedure to fix the pairing force for 
all LN and PLN calculations. However, it is well known 
that the PNP method requires another strength of the 
pairing force. Unfortunately, the average pairing gap A„ 
is not defined within the VAPNP approach, and the stan- 
dard procedure for adjusting the pairing strength is no 
longer applicable. In this study, we adjusted the VAPNP 
pairing strength to the total energy of the ^''Ca nucleus 
calculated in HFB+PLN. A much more consistent way of 
fitting the pairing strength should be based on calculat- 
ing the mass differences of the odd-mass and even-even 
nuclei, all obtained within the VAPNP method. We in- 
tend to adopt such a procedure in future applications. 

A measure of pairing correlations in a nucleus is the 
particle-particle energy (pairing energy) given by the sec- 
ond term in Eq. (|44|l . The energy of proton pairing corre- 
lations is about 2-3 MeV and it changes smoothly with N 
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NEUTRON NUMBER 

FIG. 2; Comparison between the LN, PLN, and VAPNP 
results for the chain of Ca isotopes. The upper panel shows 
the neutron pairing energies while the lower panel shows the 
total LN and PLN energies relative to the VAPNP values. 



along the isotopic chains. On the other hand, the neu- 
tron pairing is significantly affected by the shell struc- 
ture. As seen in Figs. [21 and 13 upper panels, the neu- 
tron pairing energies obtained within the LN, PLN, and 
VAPNP methods (and with pairing strengths adjusted 
as described above) are quite similar to one another. 

The lower panels of Figs. Inland |3| show differences be- 
tween the total energies obtained in the LN and PLN 
methods and those obtained in VAPNP. The LN or PLN 
results are fairly close to VAPNP for mid-shell nuclei, 
where the neutron pairing correlations are large and 
static in character. Near closed shells, pairing is dy- 
namic in nature, and the LN/PLN results deviate from 
those obtained in VAPNP. For open-shell nuclei, the PLN 
approximation is particularly good; in the calcium iso- 
topes, the deviations from the HFB+VAPNP method 
usually do not exceed 250 keV. For the closed-shell nu- 
clei, on the other hand, the LN method is not appropriate 
[T^ I20I |40| , and the energy differences increase to more 
than 1 MeV. Figures El and also show that the PLN 
method always leads to a considerable improvement over 
LN, often reducing the deviation of the total energy with 
respect to VAPNP by about 1 MeV. 
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FIG. 3: Similar to Fig.|5| except for the chain of Sn isotopes. 



As suggested in Refs. H^,!^, one can further improve 
the PLN approximation around magic nuclei by applying 
the PNP to the LN solutions obtained in the neighboring 
nuclei. This procedure is illustrated in Figs. 0] and |5l for 
the magic nuclei ^®Ca and ^^^Sn, respectively. It is seen 
that while the projection from ^^Ca nicely reproduces 
the VAPNP result in ^^Ca, the approximation fails when 
projecting from ^°Ca. Similarly, projection from the LN 
solution in ^'^"Sn (^'^^Sn) gives a better (worse) result 
than the projection of the LN solution obtained in ^^^Sn. 
We observe a similar pattern of results in other cases near 
closed shells; however, the improvement gained by pro- 
jecting from isotopes below closed shells is not sufficient 
to replace the full VAPNP calculations at closed shells. 

In order to discuss the quality of prescription to cal- 
culate the LN parameter A2 presented in Sec. IIII Bl we 
have repeated all our LN and PLN calculations with the 
effective pairing strengths G^g = aGcS scaled by factors 
of a=0.9 or 1.1 with respect to those given by Eq. (|47|l . 
In this way, we tested whether our results are sensitive to 
this phenomenological prescription. The results obtained 
for the chains of Ca and Sn isotopes are shown in Fig. 
While the LN energies (|45|l uniformly depend on the scal- 
ing factor a, the PLN energies are almost independent 
of the scaling factor. This shows that the PNP com- 
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FIG. 4: The total binding energy (with respect to a linear ref- 
erence) as a function of A*' for even-even nuclei around **Ca, 
calculated in the LN, PLN and VAPNP methods. Crosses in- 
dicate the PLN results for **Ca obtained by projecting from 
the LN solutions in neighboring nuclei ^^Ca and ^°Ca as in- 
dicated by arrows. 
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FIG. 5: Similar to Fig.QJ except for nuclei near ^^^Sn. 



ponents of the LN states weakly depend on A2 and can 
be obtained without paying too much attention to the 
way in which A2 is calculated. A rough estimate given 
by our phenomenological prescription is good enough to 
obtain reliable PLN results. On the other hand, devia- 
tions between the LN/PLN and VAPNP energies depend 
mostly on the local shell structure and visibly cannot be 
corrected by modifications of the prescription used to cal- 
culate A2. In large part, these deviations stem from the 
inapplicability of the LN/PLN method to closed-shell nu- 
clei, where the total energy in function of particle num- 
ber cannot be well approximated by the quadratic Kam- 
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lah expansion. Altogether, we conclude that the PLN 
method gives a fair approximation of the full VAPNP re- 
sults, but fails in reproducing detailed values, especially 
near closed shells. 
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FIG. 6: Total LN and PLN energies relative to the VAPNP 
values, calculated in the Ca (upper panel) and Sn (lower 
panel) isotopes with the effective pairing strengths scaled by 
a factor a. 



VI. SUMMARY AND DISCUSSION 

In this study, the variation after particle-number pro- 
jection is discussed in the context of the nuclear density 



functional theory. Specifically, we implement for the first 
time the self-consistent Skyrme HFB-I- VAPNP formal- 
ism. We demonstrate that the particle-number conserv- 
ing HFB equations with Skyrme functionals can be sim- 
ply obtained from the standard Skyrme HFB equations in 
coordinate space by replacing the intrinsic densities and 
currents by their gauge-angle-dependent counterparts. 

The calculations are carried for the Ca and Sn isotope 
chains. The VAPNP results are compared with those 
obtained with the LN and PLN methods. We demon- 
strate that the pathological behavior of LN and PLN 
methods around closed-shell nuclei can be partly cured 
by performing particle-number projection from neighbor- 
ing open-shell systems. This result is important in the 
context of large-scale microscopic mass calculations, such 
as those of ReL |22j. 

The restoration of broken symmetries in the density 
functional theory causes a number of fundamental ques- 
tions, mainly related to the density dependence of the 
underlying interaction and to the different treatment 
of particle-hole and particle-particle channels [2fil l4ll |. 
These questions and problems will be discussed in detail 
in a forthcoming paper. 



Acknowledgments 



This work was supported in part by the National 
Nuclear Security Administration under the Steward- 
ship Science Academic Alliances program through the 
U.S. Department of Energy Research Grant DE-FG03- 
03NA00083; by the U.S. Department of Energy under 
Contract Nos. DE-FG02-96ER40963 (University of Ten- 
nessee), DE-AC05-00OR22725 with UT-Battelle, LLC 
(Oak Ridge National Laboratory), DE-FG05-87ER40361 
(Joint Institute for Heavy Ion Research), DE-FG02- 
00ER41132 (Institute for Nuclear Theory); by the Polish 
Committee for Scientific Research (KBN) under Contract 
No. 1 P03B 059 27; and by the Foundation for Pohsh Sci- 
ence (FNP). 



[1] A. Bohr and B.R. Mottelson, Nuclear Structure (Ben- 
jamin, New York, 1975), Vol. II. 

[2] P. Ring and P. Schuck, The Nuclear Many-Body Problem 
( Springer- Verlag, Berlin, 1980). 

[3] D.M. Brink and R.A. Broglia, Nuclear Superfluidity: 
Pairing In Finite Systems (Cambridge Univ. Press, Cam- 
bridge, 2005). 



[4] D.J. Dean and M. Hjorth- Jensen, Rev. Mod. Phys. 75, 
607 (2003). 

[5] J. Dobaczewski and W. Nazarewicz, Phil. Trans. R. Soc. 

Loud. A 356, 2007 (1998). 
[6] J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Phys. Rev. 

108, 1175 (1957). 
[7] H. Flocard and N. Onishi, Ann. Phys. (NY) 254, 275 



11 



(1997). 

[8] A. Kamlah, Z. Phys. 216, 52 (1968). 

[9] D.C. Zheng, D.W.L. Sprung, and H. Flocard, Phys. Rev. 

C 46, 1355 (1992). 
[10] H.J. Lipkin, Ann. of Phys., 9, 272 (1960). 
[11] Y. Nogami, Phys. Rev. 134 (1964) B313. 
[12] B. GaU, P. Bonche, J. Dobaczewski, H. Flocard, and P.- 

H. Heenen, Z. Phys. A 348, 183 (1994). 
[13] P.-G. Reinhard, W. Nazarewicz, M. Bender, and J. 

Maruhn, Phys. Rev. C 53, 2776 (1996). 
[14] S. Cwiok, J. Dobaczewski, P.-H. Heenen, P. Magierski, 

and W. Nazarewicz, NucL Phys. A 611, 211 (1996). 
[15] A. Valor, J.L. Egido, and L.M. Robledo, Phys. Lett. B 

392, 249 (1997). 
[16] A. Valor, J.L. Egido, and L.M. Robledo, Nucl. Phys. A 

665, 46 (2000). 
[17] M. Bender, K. Rutz, P.-G. Reinhard, and J.A. Maruhn, 

Eur. Phys. Jour. A 8, 59 (2000). 
[18] M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz, S. Pittel, 

and D.J. Dean, Phys. Rev. C 68, 054312 (2003). 
[19] M. Anguiano, J.L. Egido, and L.M. Robledo, Phys. Lett. 

B 545, 62 (2002). 
[20] J. Dobaczewski and W. Nazarewicz, Phys. Rev. C 47, 

2418 (1993). 

[21] P. Magierski, S. Cwiok, J. Dobaczewski, and W. 

Nazarewicz, Phys. Rev. C 48, 1686 (1993). 
[22] M. Samyn, S. Goriely, M. Bender, and J.M. Pearson, 

Phys. Rev. C 70, 044309 (2004). 
[23] K. Dietrich, H.J. Mang, and J.H. Pradal, Phys. Rev. B 

135, 22 (1964). 

[24] K.W. Schmid and F. Griimmer, Rep. Progr. Phys. 50, 
731 (1987). 

[25] J.A. Sheikh and P. Ring, Nucl. Phys. A 665, 71 (2000). 
[26] M. Anguiano, J.L. Egido, and L.M. Robledo, Nucl. Phys. 

A 696, 467 (2001). 
[27] J.A. Sheikh, P. Ring, E. Lopes, and R. Rossignoli, Phys. 

Rev. C 66, 044318 (2002). 
[28] M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz, and 



P. Ring, Comput. Phys. Commun. 167, 43 (2005). 
[29] M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz, and J. 

Terasaki, Eur. Phys. J. A 25, sOl, 567 (2005). 
[30] A. Bulgac, Preprint FT-194-1980, Central Institute of 

Physics, Bucharest, 1980; nucl-th/9907088 
[31] J. Dobaczewski, H. Flocard and J. Treiner, Nucl. Phys. 

A 422, 103 (1984). 
[32] E. Perliriska, S.G. Rohozihski, J. Dobaczewski, and W. 

Nazarewicz, Phys. Rev. C 69, 014316 (2004). 
[33] Y.M. Engel, D.M. Brink, K. Goeke, S.J. Krieger, and D. 

Vautherin, Nucl. Phys. A 249, 215 (1975). 
[34] K. Hara, S. Iwasaki, and K. Tanabe, Nucl. Phys. A 332 

69 (1979). 

[35] J. Dobaczewski, W. Nazarewicz, T.R. Werner, J.-F. 
Berger, C.R. Chinn, and J. Decharge, Phys. Rev. C 53, 
2809 (1996). 

[36] J. Dobaczewski, W. Nazarewicz, and P.-G. Reinhard, 
Nucl. Phys. A 693, 361 (2001). 

[37] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and F. 
Schaeffer, Nucl. Phys. A 635, 231 (1998). 

[38] J. Dobaczewski, W. Nazarewicz, and M.V. Stoitsov, Pro- 
ceedings of the NATO Advanced Research Workshop 
The Nuclear Many-Body Problem 2001, Brijuni, Croa- 
tia, June 2-5, 2001, eds. W. Nazarewicz and D. Vretenar 
(Kluwer, Dordrecht, 2002), p. 181. 

[39] J. Dobaczewski, W. Nazarewicz, and M. V. Stoitsov, Eur. 
Phys. J. A 15, 21 (2002). 

[40] A. Valor, J.L. Egido, and L.M. Robledo, Nucl.Phys. A 
671, 189 (2000). 

[41] M. Stoitsov, J. Dobaczewski, W. Nazarewicz, P.G. Rein- 
hard, and J. Terasaki, Proc. 8th International Spring 
Seminar on Nuclear Physics, Key Topics in Nuclear 
Structure, Paestum, Italy 2004; ed. by Aldo Covello 
(World Scientific, Singapore), p. 167. 

[42] Along with the standard trace of a matrix A, 
Tr^^En in Eqs. we use Tr'A=X;„ ^nn. 



